# Load Packages -----------------------------------------------------------
library(lavaan)
library(psych)

# Load Data ---------------------------------------------------------------
source("./main/study3_cleaning.R")

# Lists to store model fit and estimated parameters
params <- list()
m_fits <- list()

# NFC ---------------------------------------------------------------------
m_nfc <- '
nfc =~ NA*nfc_teardown
+ nfc_needchaos
+ nfc_destroy
+ nfc_disastFun 
+ nfc_disastRebuild
+ nfc_burnsociety
+ nfc_burninstits
+ nfc_clearrules # rev here and down
+ nfc_upholdorder
+ nfc_disastfear # NEW
+ nfc_chaosupset # NEW
+ nfc_respectproduct # NEW
+ nfc_workinside # NEW
+ nfc_protectinstits # NEW
+ nfc_greatthings # NEW

nfc ~~ 1*nfc
'
out_s3_nfc <- cfa(m_nfc, mt)
summary(out_s3_nfc, fit.measures = T)

m_nfc_method <- '
nfc =~ NA*nfc_teardown
+ nfc_needchaos
+ nfc_destroy
+ nfc_disastFun 
+ nfc_disastRebuild
+ nfc_burnsociety
+ nfc_burninstits
+ nfc_clearrules
+ nfc_upholdorder
#+ nfc_disastfear # NEW ### comment and re-run for average loadings
+ nfc_chaosupset # NEW
+ nfc_respectproduct # NEW
+ nfc_workinside # NEW
+ nfc_protectinstits # NEW
+ nfc_greatthings # NEW

method =~ nfc_teardown
+ 1*nfc_needchaos
+ 1*nfc_destroy
+ 1*nfc_disastFun 
+ 1*nfc_disastRebuild
+ 1*nfc_burnsociety
+ 1*nfc_burninstits
+ -1*nfc_clearrules
#+ -1*nfc_disastfear # NEW ### comment and re-run for average loadings
+ -1*nfc_chaosupset # NEW
+ -1*nfc_respectproduct # NEW
+ -1*nfc_upholdorder
+ -1*nfc_workinside # NEW
+ -1*nfc_protectinstits # NEW
+ -1*nfc_greatthings # NEW

nfc ~~ 1*nfc + 0*method

# check Question Wording
nfc_needchaos ~~ nfc_destroy
'
out_s3_nfc_method <- cfa(m_nfc_method, mt)
summary(out_s3_nfc_method, fit.measures = T)
modificationindices(out_s3_nfc_method, sort = T, minimum.value = 10)

m <- out_s3_nfc_method
params[["nfc"]] <- parameterEstimates(m)[which(parameterEstimates(m)$op == "=~"),]
m_fits[["nfc"]] <- fitMeasures(m)[c("cfi", "srmr", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper")]

# Populism ----------------------------------------------------------------
m_pop <- '
pop =~ NA*pop_fewints
+ pop_crooked
+ pop_nomethink
+ pop_polsimprove # rev here down
+ pop_yesmethink
+ pop_benefitall

pop ~~ 1*pop

pop_nomethink ~~ pop_yesmethink
'
out_s3_pop <- cfa(m_pop, mt)
summary(out_s3_pop, fit.measures = T)

m_pop_method <- '
pop =~ NA*pop_fewints
+ pop_crooked
+ pop_nomethink
+ pop_polsimprove
+ pop_yesmethink
+ pop_benefitall

method =~ pop_fewints
+ 1*pop_crooked
+ 1*pop_nomethink
+ -1*pop_polsimprove
+ -1*pop_yesmethink
+ -1*pop_benefitall

pop ~~ 1*pop + 0*method

pop_nomethink ~~ pop_yesmethink
'
out_s3_pop_method <- cfa(m_pop_method, mt)
summary(out_s3_pop_method, fit.measures = T)

pars <- parameterEstimates(out_s3_pop_method)
loadings <- pars[which(pars$op == "=~" & pars$lhs != "method"), ]
mean(loadings$est[1:3])
mean(loadings$est[4:6])

m <- out_s3_pop_method
params[["pop"]] <- parameterEstimates(m)[which(parameterEstimates(m)$op == "=~"),]
m_fits[["pop"]] <- fitMeasures(m)[c("cfi", "srmr", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper")]

# Violence ----------------------------------------------------------------
m_viol <- '
viol =~ NA*viol_threatpols
+ viol_bricks
+ viol_stopbadgovt
+ viol_bullets
+ viol_noviol # rev here down
+ viol_nonviolprot # NEW
+ viol_violunaccept # NEW
+ viol_notit4tat # NEW

viol ~~ 1*viol
'
out_s3_viol <- cfa(m_viol, mt)
summary(out_s3_viol, fit.measures = T)

m_viol_method <- '
viol =~ NA*viol_threatpols
+ viol_bricks
+ viol_stopbadgovt
+ viol_bullets
+ viol_noviol
+ viol_nonviolprot # NEW
+ viol_violunaccept # NEW
+ viol_notit4tat # NEW

method =~ viol_threatpols
+ 1*viol_bricks
+ 1*viol_stopbadgovt
+ 1*viol_bullets
+ -1*viol_noviol
+ -1*viol_nonviolprot # NEW
+ -1*viol_violunaccept # NEW
+ -1*viol_notit4tat # NEW


viol ~~ 1*viol + 0*method
'
out_s3_viol_method <- cfa(m_viol_method, mt)
summary(out_s3_viol_method, fit.measures = T)

pars <- parameterEstimates(out_s3_viol_method)
loadings <- pars[which(pars$op == "=~" & pars$lhs != "method"), ]
mean(loadings$est[1:4])
mean(loadings$est[5:8])

m <- out_s3_viol_method
params[["viol"]] <- parameterEstimates(m)[which(parameterEstimates(m)$op == "=~"),]
m_fits[["viol"]] <- fitMeasures(m)[c("cfi", "srmr", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper")]

# Conspiracy --------------------------------------------------------------
m_consp <- '
consp =~ NA*consp_plots
+ consp_fewppl
+ consp_dkrun
+ consp_wars
+ consp_schoolexps # NEW, rev here down
#+ consp_fewsects # NEW
+ consp_complex # NEW
+ consp_democ # NEW
+ consp_US # NEW

consp ~~ 1*consp
'
out_s3_consp <- cfa(m_consp, mt)
summary(out_s3_consp, fit.measures = T)

m_consp_method <- '
consp =~ NA*consp_plots
+ consp_fewppl
+ consp_dkrun
+ consp_wars
+ consp_schoolexps # NEW
#+ consp_fewsects # NEW ### comment and re-run for average loadings
+ consp_complex # NEW
+ consp_democ # NEW
+ consp_US # NEW

method =~ consp_plots
+ 1*consp_fewppl
+ 1*consp_dkrun
+ 1*consp_wars
+ -1*consp_schoolexps # NEW
#+ -1*consp_fewsects # NEW ### comment and re-run for average loadings
+ -1*consp_complex # NEW
+ -1*consp_democ # NEW
+ -1*consp_US # NEW


consp ~~ 1*consp + 0*method
'
out_s3_consp_method <- cfa(m_consp_method, mt)
summary(out_s3_consp_method, fit.measures = T)

pars <- parameterEstimates(out_s3_consp_method)
loadings <- pars[which(pars$op == "=~" & pars$lhs != "method"), ]
mean(loadings$est[1:4])
mean(loadings$est[5:8])


m <- out_s3_consp_method
params[["consp"]] <- parameterEstimates(m)[which(parameterEstimates(m)$op == "=~"),]
m_fits[["consp"]] <- fitMeasures(m)[c("cfi", "srmr", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper")]

# RR ----------------------------------------------------------------------
m_rr <- '
rr =~ NA*rr_specfavr
+ rr_thard
+ rr_pdisc
+ rr_dless

rr ~~ 1*rr

rr_thard ~~ rr_specfavr
'
out_s3_rr <- cfa(m_rr, mt)
summary(out_s3_rr, fit.measures = T)

pars <- parameterEstimates(out_s3_rr)
loadings <- pars[which(pars$op == "=~" & pars$lhs != "method"), ]
mean(loadings$est[1:2])
mean(loadings$est[3:4])

m <- out_s3_rr
params[["rr"]] <- parameterEstimates(m)[which(parameterEstimates(m)$op == "=~"),]
m_fits[["rr"]] <- fitMeasures(m)[c("cfi", "srmr", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper")]


# HS ----------------------------------------------------------------------
m_hs <- '
hs =~ NA*hs_control
+ hs_exaggerate
+ hs_leash
+ hs_reasonable
+ hs_feministpower
+ hs_fewwomen

hs ~~ 1*hs

hs_reasonable ~~ hs_feministpower
'
out_s3_hs <- cfa(m_hs, mt)
summary(out_s3_hs, fit.measures = T)

m_hs_method <- '
hs =~ NA*hs_control
+ hs_exaggerate
+ hs_leash
+ hs_reasonable
+ hs_feministpower
+ hs_fewwomen

method =~ hs_control
+ 1*hs_exaggerate
+ 1*hs_leash
+ -1*hs_reasonable
+ -1*hs_feministpower
+ -1*hs_fewwomen

hs ~~ 1*hs + 0*method

hs_reasonable ~~ hs_feministpower
'
out_s3_hs_method <- cfa(m_hs_method, mt)
summary(out_s3_hs_method, fit.measures = T)

modificationIndices(out_s3_hs_method, minimum.value = 10, sort = T)


pars <- parameterEstimates(out_s3_hs_method)
loadings <- pars[which(pars$op == "=~" & pars$lhs != "method"), ]
mean(loadings$est[1:3])
mean(loadings$est[4:6])

m <- out_s3_hs_method
params[["hs"]] <- parameterEstimates(m)[which(parameterEstimates(m)$op == "=~"),]
m_fits[["hs"]] <- fitMeasures(m)[c("cfi", "srmr", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper")]


